NchooseK_FlatBranchHybrid <- function(BD,K,BranchProbs,CompProbs) {
#   def subset_sum(numbers,target):
  
  # See Wikipedia for description
  #
  # Inputs
  #   N - scalar - number of bins
  #   BD - 1xN vector - gives the bin depths, thereby exclusing certain combinations
  #   K - scalar - number of objects fitting into those bins
  #   BranchProbs - 1xN list w/ asymmetrical vectors - 1-K probabilities of each branch
  #   CompProbs - 1xN vector - Probabilities of the parent components of those branches
  #
  # Outputs
  #   Combos - List vector - provides all the combos
  #
  # Written by: Sean M. Gonzalez, Willcor Inc., December 2012
  #______________________________________________________________
  
  pcird <- file("ProbCompIsReplaced.data", "a")
  
  # we need an intermediate function to start the recursion.
  # the recursion start with an empty list as partial solution.
  if (any(BD==0)) { stop("A bin can't have zero depth!") }
  if (K>sum(BD)) { cat("There is no possible solution for this K =",K,"in these bins =",BD,"\n",file=pcird) }
  
  # I need to begin with a set of number combinations that add to K
#   BD = sort(BD) # I can't sort because then it's out of index with it's corresponding list vector
  N = length(BD)
  numbers = 1:K
  nx = numbers > max(BD)
  numbers = numbers[!nx]
  target = K
  PCNcK = list("Prob"=0,"NumEvents"=0)
  cat("\nInitial PCNcK:",file=pcird)
  cat("\nProb =",     PCNcK$Prob,     file=pcird)
  cat("\nNumEvents =",PCNcK$NumEvents,"\n",file=pcird)
  partial = matrix(nrow=0,ncol=0)
  
  PCNcK <- subset_sum_recursive(N,BD,numbers,target,partial,BranchProbs,CompProbs,PCNcK)
  
  cat("\nLeaving Stars & Bars\n",file=pcird)
  return(PCNcK)
}  
#==================================================================
subset_sum_recursive <- function(N,BD,numbers,target,partial,BranchProbs,CompProbs,PCNcK) {
#   http://stackoverflow.com/questions/4632322/finding-all-possible-combinations-of-numbers-to-reach-a-given-sum
  pcird <- file("ProbCompIsReplaced.data", "a")
  sumPartial = sum(partial)
  lBD = length(BD) # Wish I could have a persistant variable...
  
  # Check if the partial sum is equals to target
  if (sumPartial == target) { 
    lP  = length(partial)
    if (lP !=lBD) {
      delta = lBD - lP
      if (delta <0) { 
        cat("partial =",partial,": BD =",BD,"\n",file=pcird)
        stop("How is number set greater than the bins available?")}
      partial = c(partial,integer(delta))
    } # Pad partial with zeroes
    
    cat("We've reached our target! \n",file=pcird)
    
    PCNcK <- CountBinPermutations(partial,BD,BranchProbs,CompProbs,PCNcK)
    
    cat("PCNcK w/ partial =",partial,"\n",file=pcird)
    cat("PCNcK$NumEvents =",PCNcK$NumEvents,"\n",file=pcird)
    cat("PCNcK$Prob =",     PCNcK$Prob,     "\n",file=pcird)
    return(PCNcK)
  }
  if ((sumPartial > target) || length(partial)==lBD){
    cat("Sum was greater than target!?\n",file=pcird)
    cat("Partial =",partial,": target =",target,"\n",file=pcird)
    cat("Partial =",partial,": BD =",BD,"\n",file=pcird)
    return(PCNcK) # if we reach the number why bother to continue
  }
  
  cat("\nLooking to drop some levels.\n",file=pcird)
  cat("N =",N,": numbers available =",numbers,"\n",file=pcird)
  if ((N >0) && (length(numbers)>0)) {
    
    for (i in 1:length(numbers)) { 
      # We're adding this number to what we're given, 
      # and removing it from the list for the downstream recursions
      n = numbers[i]
      remNums = numbers[-i] # Remaining Numbers to play with
      cat("\nEvaluating this number to add:",n,"\n",file=pcird)
      
      # First: Determine how many bins this first number occupies: 
      # binsAvail, binsPossible
      {
      binsPossible = which(BD >=n)
      numBinsPoss = length(binsPossible)
      numBinsReq = floor(target/n)
      cat("BD =",BD,": numbers =",numbers,": n =",n,": target =",target,": numBinsPoss =",numBinsPoss,": numBinsReq =",numBinsReq,"\n",file=pcird)
      if (numBinsPoss ==0) { break } # numbers is sorted, so if not this then none
      if (numBinsReq >numBinsPoss) {
        numBinsReq = numBinsPoss # We'll make up the remainder with another number
      } }
      # Can we succeed with what we're given: baIdx
      {
      baIdx = 1:lBD # bins available Idx
      baIdx = baIdx[-(lBD:(lBD-numBinsPoss))] # [-(1:numBinsPoss)] # BD is sorted, remember
      if (N*max(numbers) < target) {
        return(PCNcK) # Can't succeed with what we're given
      } }
      
      # Use the same number as many times as possible 
      #   given what's remaining at this recursion level
      newPartial = partial
      for (bo in 1:numBinsPoss) {
        cat("\nUsing the number ",n," for the ",bo," time\n",file=pcird)
        
        thisTotal = bo*n + sumPartial
        newPartial = sort(t(as.matrix(c(newPartial,n))))
        binsLeft = (N - bo)
        if (length(newPartial)>lBD) {
          cat("Trying to use too many number to reach the goal.\n",file=pcird)
          cat("newPartial =",newPartial,": BD =",BD,"\n",file=pcird)
          return(PCNcK)
        }
        { # I need to make sure that it's feasible to add this number given numbers available.
          sBD = sort(BD)
          #     srebmun = sort(numbers,decreasing=TRUE)
          sNumbers = sort(newPartial)    
          currentDepth = 1 # The numbers reach the first BD
          cat("sBD = ",sBD,": sNumbers = ",sNumbers,"\n",file=pcird)
          for (m in 1:length(sNumbers)) {
            num = sNumbers[m]
            if ((num ==0) || (num ==1)) { next }
            thisDepth = min(which(sBD >= sNumbers[m])) -m+1
            cat("thisDepth =",thisDepth,"\n",file=pcird)
            currentDepth = max(currentDepth, thisDepth)
          }
          if (currentDepth >=lBD) {
            return(PCNcK) # You can't add anything new, and the rest are even bigger.
          } } # Can we add anything new?

        cat("\nGoing down another level! \n",file=pcird)
        cat("Partial vs. newPartial =",partial,"vs.",newPartial,"\n",file=pcird)
        cat("PCNcK$NumEvents =",PCNcK$NumEvents,"\n",file=pcird)
        cat("PCNcK$Prob =",     PCNcK$Prob,     "\n",file=pcird)
        
        PCNcK <- subset_sum_recursive(binsLeft,BD,remNums,target,newPartial,BranchProbs,CompProbs,PCNcK) 
        
        cat("\nBacked out a level! \n",file=pcird)
        cat("Partial vs. newPartial =",partial,"vs.",newPartial,"\n",file=pcird)
        cat("PCNcK$NumEvents =",PCNcK$NumEvents,"\n",file=pcird)
        cat("PCNcK$Prob =",     PCNcK$Prob,     "\n",file=pcird)
      }
      cat("Used Number ",n," in ",numBinsPoss," possible bins.\n",file=pcird)
    }
  }
  else {
    # You have no more options with this permutation
  }
  return(PCNcK)
}
#===============================================================================
CountBinPermutations <- function(numbers,BD,BranchProbs,CompProbs,PCNcK) {
  
  # I don't need to assert the values of the inputs, such as them adding to the target
  #   or the length of numbers being equal to BD, because 
  #   that's asserted before calling this function
  
  pcird <- file("ProbCompIsReplaced.data", "a")
  cat("\nInside 'CountBinPermutations'\n",file=pcird)
  
  # Initialize
  lBD = length(BD)
  uniqueNums = unique(numbers)
  lUN = length(uniqueNums)
  lNZN = sum(numbers>0) # Length of Non-Zero-Numbers
  binsTaken = 0
  binPerms  = 1
  binsSoFar = 0
  failPerms = numeric(lUN-1)
  failProbs = failPerms
  passProbs = 1
  
  count = 0
  for (i in 1:lUN) {
    n = uniqueNums[i]
    # bx is the bin indices
    if (n==0) {
      bx = 1:lBD } #!logical(lBD) } # Assign bx
    else {
      bx = which(BD >=n)
    }
    nx = which(numbers==n)
    potBins    = length(bx)
    binsNeeded = length(nx)
    # I can use bx & nx to get a more precise cost average
    
    if ((n!=0) && (binsNeeded > potBins)) { 
      cat("numbers =",numbers,": BD =",BD,": this# =",n,"\n",file=pcird)
      stop("Should not have been able to get this combination of numbers passed in...") }
    if ((n==0) && (potBins !=length(BD))) {
      cat("potBins =",potBins,": BD =",BD,": this# =",n,"\n",file=pcird)
      stop("Zero fits everywhere...")
    }
    
    binsTaken = binsTaken + potBins
#     PartFailProb[k] = sum(combn(C,k,function(kx) {xx[kx]=TRUE; prod(ChildProbs[xx,k])*prod((1-ChildProbs[!xx,k]))}))
    
    # Get the average probability of the bins this number may occupy
    # The one exception is 0 (zero), 
    #   which means that branch survives, and uses the Parent probability
    if (n==0) { 
      if (potBins !=length(BD)) { 
        stop("If there needs to be a passing component, then any component could pass") }
      # Don't need number of ways remaining components can pass, it's always one!!
      passProbs = mean(1-CompProbs) } # All input probs are failing probs.
    else { 
      count = count +1
      failPerms[count] = choose(potBins,binsNeeded)
      cat("BD =",BD,": bin indices (bx) =",bx,": n=",n,"\n",file=pcird)
      cat("\nfailProbs to average:\n",file=pcird)
      
      temp = numeric(length(bx))
      for (bb in 1:length(bx)) {
        xx = bx[bb] # bin index
        cat("temp =",temp,": mean(temp) =",mean(temp),": Branch Length =",length(BranchProbs[[xx]]),"\n",file=pcird)
        cat("UniqueNum =",n,": bin index =",xx,": bin indices =",bx,": bin step =",bb,"\n",file=pcird)
        cat("BranchProbs =",BranchProbs[[xx]],": Number of Branches =",length(BranchProbs),"\n",file=pcird)
        temp[bb] = BranchProbs[[xx]][n] # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#       cat("BranchProbs =",BranchProbs,"\n") # !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      }
      cat("temp =",temp,"\n",file=pcird)
      failProbs[count] = mean(temp)
    }
  }
  numEvents = prod(failPerms)
  totFailProb = (numEvents *mean(failProbs)^lNZN *mean(passProbs)^(lBD-lNZN))
  cat("numbers =",numbers,"\n",file=pcird)
  cat("BD =",BD,"\n",file=pcird)
  cat("numEvents =",numEvents,"\n",file=pcird)
  cat("failPerms =",failPerms,"\n",file=pcird)
  cat("failProbs =",failProbs,"\n",file=pcird)
  cat("passProbs =",passProbs,"\n",file=pcird)
  PCNcK$NumEvents = PCNcK$NumEvents + numEvents
  PCNcK$Prob      = PCNcK$Prob      + totFailProb
  cat("PCNcK$NumEvents =",PCNcK$NumEvents,"\n",file=pcird)
  cat("PCNcK$Prob =",     PCNcK$Prob,": This Prob =",totFailProb,"\n",file=pcird)
  return(PCNcK)
}